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Abstract 

The  description  of  the  total  energy  of  a  solid  in  terms  of  the  effective  medium  theory 
has  been  studied.  Using  the  density  functional  theory  an  exact  symmetric  formula  is 
presented  for  describing  the  total  energy  in  terms  of  embedding  energies  of  individual 
atoms.  Different  approximate  formulae  are  derived  and  compared  to  the  conventional 
pair  potential  representation,  to  the  model  of  Gordon  and  Kim,  and  to  the  semiempiri- 
cal  embedded  atom  scheme  of  Daw  and  Baskes.  The  theory  is  applied  to  estimate  the 
pair  potential  in  rare  gases  He,  Ne,  and  Ar,  and  the  cohesive  properties  of  A1  metal.  A 
fair  agreement  with  the  experimental  results  is  obtained.  fj 
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1.  Introduction 


The  study  of  static  and  dynamic  properties  of  lattice  defects  using  molecular  dynamics 
is  based  on  the  pair  potential  representation  of  the  structural  part  of  the  total  energy  of 
the  system1.  In  metals  this  is  usually  written  as 

Etot  =  yEXW1*,)  +  Evob  (1-1) 

^  i  j*i 

where  Vp  is  the  pair  potential,  Rf  a  position  of  an  atomic  nucleus  and  N  is  the  total 
number  of  atoms.  Evol  is  an  energy  term  which  depends  only  on  the  volume  (or  den¬ 
sity)  of  the  system  but  not  on  the  structure.  A  pair  potential  expression  of  the  form  of 
Eq.  (1.1)  is  formally  obtained  by  describing  the  metal  ions  by  pseudopotentials  and 
applying  second  order  perturbation  theory  to  calculate  the  total  energy2-3.  Several 
different  schemes  have  been  proposed  to  determine  the  pseudopotential  and  the 
corresponding  pair  potential.  In  alkali  and  other  nontransition  metals  the  results 
obtained  for  example  for  the  phonon  spectra  are  in  fair  agreement  with  the  experimen¬ 
tal  results4.  However,  the  pair  potential  picture  becomes  questionable  when  applied  to 
lattiCe  defects  which  cause  large  perturbations  to  the  nearly  homogeneous  valence  elec¬ 
tron  density  of  the  perfect  metal.  Also  the  treatment  if  impurities  in  the  pseudopoten¬ 
tial  scheme  is  not  straightforward. 

Recently,  another  method  has  been  proposed  for  calculating  the  interionic  interac¬ 
tions5,6.  In  this  so-called  "embedded  atom"  method  one  tries  to  write  the  total  energy 
in  the  form 

Etot  =  E  Vi))  +  tI  WrR,)>  (1.2) 

i  L  i*j 

where  F  is  a  function  of  the  electron  density  ni  at  the  site  of  the  atom  i  and  O  is  the 
residual  part  of  the  pair  ion  interaction  which  can  not  be  included  in  the  first  term  of 
Eq.  (1.2).  The  density  «,  is  the  electron  density  at  site  ft,  when  the  atom  from  that  site 
is  removed,  i.e.  it  is  the  electron  density  provided  by  the  surrounding  atoms.  The  idea 
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is  to  approximate  the  total  electron  density  as  a  superposition  of  atomic  densities  na  in 
which  case  ni  is 

N 

«,(R«)  =  Z  na(Rr^j)-  (1.3) 

j*‘ 

The  equations  (1.2)  and  (1.3)  are  essentially  as  easy  to  use  in  the  computer  simulations 
as  the  standard  pair  potential,  Eq.  (1.1). 

The  function  F(n )  is  generally  nonlinear  and  then  the  first  term  implicitly  takes  into 
account  the  volume  dependent  terms  of  Eq.  (1.1).  Also,  the  expression  (1.2)  includes 
multi-ion  interactions18.  It  is  expected  that  the  embedded  atom  scheme  can  then  better 
describe  defects  like  vacancies  and  surfaces  where  the  density  variations  are  large. 

The  idea  behind  the  Eq.  (1.2)  comes  from  the  so-called  "effective  medium"  or  "pseu¬ 
doatom"  theory  first  presented  for  calculating  the  energetics  of  single  impurities  in 
metals7,8.  Daw  and  Baskes5  have  applied  successfully  the  embedded  atom  method  for 
many  transition  metal  systems  by  taking  an  empirical  approach,  the  functions  F  and  O 
are  fitted  to  reproduce  some  experimental  data. 

The  purpose  of  the  present  paper  is  to  look  more  formally  the  possibilities  of  deriving 
interionic  interactions  of  the  form  of  Eq.  (1.2)  using  the  effective  medium  theory. 
First  the  total  energy  of  the  solid  (in  any  ionic  arrangement)  is  written  in  a  symmetric 
way  in  terms  of  the  electron  densities  provided  by  the  surrounding  atoms.  This  can  be 
done  formally  exactly  using  the  density  functional  theory.  Then  approximate  formulae 
are  derived  by  using  the  same  approaches  which  have  been  earlier  applied  for  single 
impurities.  It  is  shown  that  the  first  order  approximation,  in  the  limit  of  a  weak  distur¬ 
bance  from  a  homogeneous  lattice,  reduces  to  the  conventional  pair  potential  picture 
with  the  same  pair  potential  as  obtained  from  the  perturbation  theory.  In  the  case  of 
rare  gases,  on  the  other  hand,  the  result  is  closely  related  to  the  Gorgon-Kim9  model 
for  molecular  binding. 


In  Section  2  we  present  a  formal  derivation  of  the  total  energy  in  terms  of  density 
functionals.  In  Section  3  approximations  are  made  for  obtaining  a  usable  equation  of 
the  form  of  Eq.  (1.2).  In  Section  4  the  results  are  compared  to  the  conventional  pair 
potential  picture,  to  the  scheme  of  Gordon  and  Kim9,  and  to  the  semi-empirical  for¬ 
malism  of  Daw  and  Baskes5.  Approximative  (  ab  initio  )  results  for  rare  gases  and  A1 
metal  are  given  in  Section  5.  Section  6  includes  discussion  and  conclusions. 


-  V- 


2.  Formal  theory 

According  to  the  density  functional  theory10,11  the  total  energy  of  system  of  electrons 
in  an  external  potential  can  be  expressed  as  a  functional  of  the  total  electron  density. 
This  functional  has  the  variational  property  that  the  exact  ground  state  electron  density 
minimizes  it12.  In  a  solid  the  external  potential  is  provided  by  the  atomic  nuclei  which 
are  fixed  a  sites  /?,•  (adiabatic  Bom-Oppenheimer  approximation  is  assumed).  The  total 
energy  of  the  metal  consisting  of  N  atoms  is 

E%T  =  ENR[n)  (2.1) 

where  n  is  the  total  ground  state  electron  density  and  the  subscript  R  indicates  that  the 
functional  depends  on  the  sites  R,  and  charges  Z,  of  the  nuclei  (also  the  electrostatic 
nucleus-nucleus  repulsion  is  included).  If  one  atom  is  removed  from  the  metal  the 
energy  change,  the  so-called  embedding  energy,  is 

AEi  =  Ej0T  -  Ejo\  -  Eatom,  (2.2) 

where  EATO M  is  the  total  energy  of  a  free  atom  in  vacuum.  The  subscript  i  indicates 
that  the  energy  depends  upon  which  atom  is  removed.  The  embedding  energy  can  be 
written  as7 

A£,  =  A £,[«,],  (2.3) 

where  n-t  is  the  self-consistent  ground  state  electron  density  in  the  system  after  the 
atom  i  has  been  removed  (but  other  atom  sites  are  kept  fixed).  Equation  (2.3)  is  a 
trivial  generalization  of  the  density  functional  theorem.  The  density  «;  uniquely 
defines  the  external  potential  of  the  N- 1  nuclei.  Then,  since  the  site  of  the  removed 
atom  R,  is  known,  also  the  external  potential  of  the  original  system  with  N  atoms  is 
known  and  both  Ejq\ ■  and  E% T  can  be  formally  expressed  as  functionals  of  the  same 
density  This  is  the  key  idea  of  the  effective  medium  theory  when  applied  to  calcu¬ 
late  energetics  of  impurities  in  metals.  It  simply  means  that  the  binding  energy  of  the 


C  I 


a 

i 


5& 

•Sia 

.‘M 


impurity  can  be  calculated  directly  from  the  unperturbed,  electron  density  of  the  host 
metal7,8. 

Using  Equations  (2.2)  and  (2.3)  we  write  the  total  energy  of  the  metal  now  in  terms  of 
the  embedding  energies.  By  removing  one  atom  we  can  immediately  write 

ETot  =  +  +  Eatom.  (2.4) 

We  can  successively  remove  atoms  one  by  one  and  get  finally 

Etot  =  ME  atom  +  S  A£if/l123...i]  (2-5) 

1=1 

where  N  is  the  number  of  atoms  in  the  system  and  nl2j  j  is  the  self-consistent  electron 
density  of  the  system  where  atoms  at  the  sites  1,2,3,  •  •  •  /  have  been  removed.  Each 
term  of  this  sum  depends  upon  which  atom  is  removed  and  which  atoms  had  been 
taken  out  before.  By  taking  an  average  of  all  possible  orders  of  removing  the  atoms 
the  equation  (2.6)  can  be  made  symmetric: 

Etot  =  NEAT0M  +  jj 

N(N-[)(N-2)  k  ljU 

=  NEAtom  +  +  S2  +  S3  +  •  •  ■  (2.6) 

The  prime  in  the  sums  above  indicates  that  terms  which  have  two  or  more  indexes  the 
same  are  omitted.  Each  sum  in  Eq.  (2.6)  contributes  about  the  same  order  of  magni¬ 
tude  to  the  total  energy.  To  find  a  convergent  series  we  first  define  sums 
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In  these  sums  only  those  terms  in  which  each  of  the  removed  atoms  are  close  to  each 
other  are  non-negligible.  If  the  atoms  are  far  apart  they  do  not  interact 
(A Ej[nij]  =  AEj[nj]  if  |R— Rj\  large).  The  sums  S  are  related  to  the  sums  S  in  Eq.  (2.6) 


5n  “  ^  +  L-2)^-1 


+  •  •  •  +5,. 


Equation  (2.6)  can  now  be  rewritten  as 
ETOt  =  neatom  + 


+  yZ:.(a Ej{nij\  -  A£,[n,]] 

+  jZ.'  Ja£,[^]  -  A Ek[nik\  -  A Ek[njk]  +  A £*[«*]] 


(2.10) 


Now  in  each  sum  (except  in  the  first)  each  term  is  negligible  if  any  pair  of  the  sites 
are  far  from  each  other.  It  is  then  expected  that  the  expansion  converges  fast  even  if 
each  sum  has  seemingly  more  terms  than  the  previous  one.  (For  example  later  it  will 
be  shown  that  in  calculating  electrostatic  interactions  between  atoms,  only  the  first  two 
sums  in  this  expansion  are  nonzero). 

In  this  functional  form  the  Eq.  (2.10)  is  exact.  In  a  large  perfect  crystal,  each  atom  is 
in  a  similar  surrounding.  The  cohesive  energy  can  the  be  expressed  as 


=  -A£,K1  -  -  A£,H)] 


(2.11) 


and  since,  by  definition,  A £,[«,]  is  the  embedding  energy  of  an  atom  in  a  pre-existing 
vacancy,  the  vacancy  formation  energy  (for  an  unrelaxed  vacancy)  is 


Evac  =  "  A£,T«yl] 


(2.12) 
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The  exact  symmetric  expression,  Eq.  (2.10),  is  now  a  starting  point  in  trying  to  find 

applicable  approximate  formulae.  I 


3.  Approximative  expressions 

3.1.  Electron  density 

The  electron  density  which  determines  the  embedding  energy  in  the  functional 
A£,[«i23.j]  is  the  self-consistent  ground  state  electron  density  of  the  system  where 
atoms  from  sites  123. ..i  are  removed.  It  is  obvious  that  any  expression  based  on  the 
embedding  energy  functionals  can  not  be  useful  if  the  self-consistent  electron  density 
has  to  be  calculated  for  each  atomic  configuration.  (It  would  then  be  easier  to  calcu¬ 
late  directly  the  total  energy  of  the  whole  system  self-consistently  and  the  effective 
medium  theory  would  not  be  needed  at  all).  For  making  a  practical  scheme,  the  elec¬ 
tron  density  has  to  be  approximated  by  a  superposition  of  densities  associated  with 
each  atom.  These  densities  /ia(r-R;)  do  not  necessarily  have  to  be  the  densities  of  free 
atoms  but  can  be  densities  of  some  kind  of  pseudoatoms  which  better  describe  the 
screening  charge  in  the  metal.  (In  the  linear  screening  of  pseudopotentials  na  would  be 
exactly  the  screening  charge  of  the  pseudopotential). 

In  the  following  we  will  always  make  the  approximation  that  the  density  can  be 
expressed  as 

N 

«123.j(r)  =  E  «a(r-R j),  (3.1) 

y=‘+i 

and  specify  later,  whenever  necessary,  what  exactly  is  meant  by  the  pseudoatom  den¬ 
sity  na. 

3.2.  Local  approximation  for  the  embedding  functional 

As  an  illustrative  example  we  will  now  assume  that  the  embedding  energy  is  a  func¬ 
tion  of  the  local  electron  density.  This  is  the  simplest  form  of  the  embedding  energy 
functionals.  The  function  can  be  calculated  for  example  by  embedding  the  atom  in  a 
homogeneous  electron  gas,  in  which  case  it  can  be  viewed  as  the  first  term  in  a 
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systematic  expansion  of  the  embedding  energy  in  terms  of  density  gradients  or  pertur¬ 
bation  corrections7.  In  the  local  approximation 

AE.-K1  =  A£*om(nI<Ri)),  (3-2) 


where  A Ehom(n)  is  the  embedding  energy  of  the  atom  in  a  homogeneous  electron  gas 
of  density  n.  To  simplify  the  notations  we  have  dropped  the  subscript  i  from  the  func¬ 
tion  A Ehom  assuming  that  there  is  only  one  kind  of  atoms  in  the  system.  Using  Eq. 
(3.1)  for  the  density  we  can  write 

«y(R,)  =  «/R;)-«a(RrR<)  (3-3) 

and  the  function  A Ehom(ni})  needed  in  Eq.  (2.10)  can  be  expanded  as 

.  .  3A Ehom(n:) 

A Ehom{n  j  =  A Ehom{nj)  -  na{ RrR;) - — +  •  •  •  (3.4) 

Expanding  also  AEhom(nijk)  in  the  same  way  and  substituting  in  Eq.  (2.10)  we  find  (by 
neglecting  a  small  second  order  term  and  all  higher  order  terms,  see  Appendix  A) 


Ejot  -  NEatom  +  Yi 


,  i  dAE^n:)  i  d2AEhom(ni) 

^  ("f) '  + 


(3.5) 


By  defining  a  function 


cAom.  .  _  p  .  _  1  +  i„2  3l4£'"”M 

1  ^n)~tATOM  +  At  (")-—«  JZ 


dn2 


(3.6) 


the  total  energy  of  the  metal  can  be  written  as 

£tot  =  Z^K-fR.-))-  (3-7) 

i=i 

This  is  of  the  desired  form  of  Eq.  (1.2).  In  this  simplest  local  approximation  there  is 
no  additional  pair  potential  sum.  The  relationship  (3.6)  between  the  function  Fhom  and 
the  embedding  energy  function  A Ehom  can  also  be  obtained  directly  by  requiring  that 
the  total  energy  is  written  in  the  form  of  Eq.  (2.7)  as  shown  in  the  Appendix  A.  A 
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purely  local  expression  of  the  form  of  Eq.  (3.7)  can  not  be  a  good  approximation  for 
solids  since  it  results  unreasonable  results  for  elastic  constants  as  shown  by  Daw  and 
Baskes5. 
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3.3.  Corrections  from  the  perturbation  theory 

One  systematic  way  to  improve  the  local  approximation  in  the  effective  medium  theory 
is  to  use  the  perturbation  theory7,8.  The  first  order  correction  to  the  homogeneous  elec¬ 
tron  gas  term  of  Eq.  (3.2)  is 

AE^fR.)  =  J^rApfr-R^Sv^V)  (3.8) 

where  Ap,(r)  is  the  atom  induces  charge  density  in  the  homogeneous  electron  gas 
(Ap(r)  =  An(r)-ZS(r),  where  An  is  the  induced  density)  and  8vexf(r)  is  the  difference  in 
the  external  potentials  between  the  homogeneous  electron  gas  and  the  real  metal. 
Equation  (3.8)  can  be  written  in  a  more  useful  form  by  using  for  the  homogeneous 
electron  gas  density  an  average  density 

1  .  ,  ,  ,  Ap(r'-R.) 

ni<Ri)  =  -~\d  rni(r)jd3r'-— r--~  (3.9) 

where  a  is  defined  as 


(3.10) 


and  is  assumed  to  be  nonzero  (which  is  not  necessarily  true  for  all  Ap7,8).  Using  Eqs. 
(3.9)  and  (3.10)  it  is  straightforward  to  show  that  the  first  order  correction  reduces  to 
the  electrostatic  interaction  between  the  induced  charge  density  of  the  embedded  atom, 
Ap,  and  the  total  charge  density  of  the  system  in  which  the  atom  is  embedded.  Eq. 


(3.8)  then  becomes 


AE<1>(RI)  =  JAAp(r-R,)<Xr) 


(3.11) 


where  0  is  the  total  electrostatic  potential  of  the  system  (without  the  potential  of  the 
atom  to  be  embedded).  To  be  consistent  with  the  approximation  (3.1)  the  pseudoatom 
density  na  should  be  taken  to  be  the  electron  density  corresponding  to  the  induced 
charge  Ap  (i.e.  Ap  minus  the  nuclear  charge).  With  this  approximation  the  correction 
A£*^  can  be  written  as  follows 


...  ,  ,  .  Ap(r-R;)Ap(r'-R.) 

A£«>(R <)  =  Z  jcl'rfar'  ^ 


(3.12) 


Substituting  this  to  Eq.  (2.10)  it  turns  out  that  only  the  two  first  sums  of  the  expression 
are  non-zero.  Including  also  the  local  part  of  the  embedding  energy  from  the  previous 
section  the  total  energy  of  the  metal  will  be 


Etot  =  £f‘“"«Ri))  +  ±TlS<Prj<Pi' 


I 


,  Ap(r-R,)Ap(r-R7) 
ir-r'| 


(3.13) 


This  expression  is  again  of  the  form  of  Eq.  (1.2),  the  pair  potential  being  now  the  elec¬ 
trostatic  interaction  between  the  pseudoatoms.  The  fact  that  the  average  density  n  now 
appears  in  the  function  Fhom  does  not  make  the  application  of  the  Eq.  (3.13)  more 
difficult  since  the  averaging  can  be  made  to  the  pseudoatom  density  before  the  summa¬ 
tion  of  the  total  density  (i.e.  the  pseudoatom  density  in  Eq.  (3.1)  is  replaced  by  an 
averaged  pseudoatom  density  na  ).  However,  there  is  a  self-consistency  requirement: 
the  induced  density  Ap  should  be  calculated  by  embedding  the  atom  in  a  homogeneous 
electron  gas  of  density  n  which  depends  on  Ap  through  Eqs.  (3.1)  and  (3.9).  This  con¬ 
dition  can  not  be  strictly  fulfilled  if  the  superposition  approximation  for  the  total  elec¬ 
tron  density  is  required.  The  most  consistent  way  is  to  calculate  Ap  in  a  density  TTq 
which  is  an  average  over  all  atom  sites  (as  defined  in  Appendix  C,  Eq.  (C.3)).  The 
induced  charge  density  is  spherically  symmetric.  At  small  distances  from  the  nucleus  it 
is  governed  by  the  core  electrons  and  is  very  close  to  that  of  a  free  atom.  At  large 
distances  Ap(r)  has  Friedel  oscillations  whereas  the  free  atom  density  goes  to  zero 
exponentially.  If  the  short  range  atomic  like  behavior  is  dominating  in  the  electrostatic 
sum  of  Eq.  (3.13)  then  the  results  should  not  be  sensitive  on  which  density  Ap  is 
determined.  However,  if  the  long  range  Friedel  oscillations  turn  out  to  be  important, 
then  the  self-consistency  becomes  essential  and  also  the  approximation  of  the  superpo¬ 
sition  of  the  pseudoatom  densities  questionable. 


In  the  perturbation  theory  also  higher  order  corrections  can  be  systematically  added  in 
calculating  the  embedding  energy  functional  in  a  inhomogeneous  electron  gas.  Since 
we  want  to  keep  the  approximation  of  Eq.  (3.1)  there  is  not  any  obvious  systematic 
way  to  improve  the  total  energy  functional  ,  Eq.  (2.10),  using  higher  order  perturbation 


3.4.  Corrections  using  the  mixed  perturbation  scheme 


For  taking  advantage  of  the  local  nature  of  the  screening  in  metals  and  for  avoiding 
the  problems  arising  from  the  long  range  Friedel  oscillations,  Ndrskov  13  has  proposed 
a  "mixed  perturbation  scheme".  The  key  idea  is  that  in  the  close  vicinity  of  the  embed¬ 
ded  atom  (region  a,)  the  potential  is  governed  by  the  atom  and  can  be  approximated  to 
be  the  same  both  in  the  homogeneous  electron  gas  and  in  the  real  metal.  Further  out, 
on  the  other  hand,  the  potential  can  be  approximated  to  be  that  of  the  host,  i.e. 
unaffected  by  the  embedded  atom.  The  lowest  order  correction  to  AEhom  in  this 
approach  can  be  written  as  follows6,13, 

8A Ei  =  Jd3rAp(r— R,)8<jra,(r) 

Oi 

+  JdVApCr-R^Sv^r) 

a, 

Ef 

+  5 j  deeAn(e)  (3.14) 


where  S<Jfa,’(r)  is  the  part  of  the  electrostatic  potential  of  the  metal  which  is  caused  by 
the  charges  outside  the  region  a4-  and  5v^,(r)  is  the  change  in  the  external  potential 
caused  by  the  positive  charges  inside  the  region  a,-.  The  last  integral  comes  from  the 
change  in  the  one-electron  energy  eigenvalues  when  the  atom  is  moved  from  the 
homogeneous  electron  gas  to  the  real  metal.  The  first  two  integrals  are  extended  over 
a  sphere  at-  centered  ar  R(.  In  practice  the  sphere  a,-  is  always  so  small  that  inside  it 
there  is  no  other  nuclei  than  that  of  the  embedded  atom.  Then  Sv^,  is  caused  by  the 
positive  background  charge  of  the  homogeneous  electron  gas.  By  making  again  the 
approximation  (3.1). of  superposition  of  pseudoatom  densities,  Eq.  (3.14)  can  be  written 
as  (see  Appendix  B  for  details) 


5A£(  =  £jd3rj</V 


Ap(r-R,)Ap(r-R;) 

Ir-r'l 


t.  r  ,,  rj,  ,Ap(r-R,)Ap(r-R) 

-X  j  <fr\db - — - 

j*i  -a,  11  1 


5  J  deeAn(e) 


(3.15) 


where  the  notation  - a ,•  under  the  integral  means  that  the  sphere  omitted  in  the 
integration. 

The  change  in  the  one-electron  eigenvalues  is  difficult  to  estimate.  Ndrskov  et  al6 
have  studied  approximate  ways  to  derive  this  term  and  thus  go  beyond  the  result  of 
Eq.  (3.13).  Here  we  only  want  to  mention  that  if  the  potential  outside  the  sphere  a,-  is 
assumed  to  be  purely  electrostatic  and  apply  simple  perturbation  theory  to  estimate  the 
change  in  the  eigenvalues,  the  last  term  in  Eq.  (3.15)  cancels  exactly  the  second  term. 
The  result  of  the  normal  perturbation  theory  of  Section  3.3  is  then  recovered.  Note  that 
in  this  limit  the  result  is  independent  of  the  choice  of  the  radius  of  the  sphere  a-.. 


4.  Relation  to  other  methods 


4.1.  Pair  potential 


The  present  result  of  the  effective  medium  theory,  Eq.  (3.13),  which  includes  perturba¬ 
tion  correction  to  the  local  energy  functional,  can  be  related  to  the  conventional  pair 
potential  picture  derived  from  the  pseudopotentials.  This  can  be  done  by  extracting 
from  the  first  term  of  Eq.  (3.13)  the  pair  interaction  part  and  the  structure  independent 
part.  The  straightforward  derivation  is  done  in  Appendix  C  and  it  gives 


ETOT  ~  sj EATOM  +  A Ehom(n0)  -  COi0j 


[Ap(r-RI)-na(r-R1)]Ap(r,-R/) 

T  IlP'P' - - 


(4.1) 


i  j*i 


where  nQ  is  an  average  of  the  (averaged)  densities  in  Eq.  (3.9)  and  Ap  and  na  are  the 
charge  density  and  electron  density  of  the  pseudoatom,  respectively.  The  first  sum  is 
now  structure  independent  but  depends  of  the  volume  of  the  system  through  n0  .  The 
second  sum  is  the  pair  interaction:  it  is  the  sum  of  the  electrostatic  interactions 
between  a  screened  atom  (screened  pseudopotential)  and  a  bare  pseudopotential 
expressed  in  the  form 

-  r  .3  ,  AP(r-Q  ~  ”a(r-r/) 


y  pseudo 


-K 


|r-r'| 


(4.2) 


Note  that  na  is  here  the  electron  density  corresponding  to  the  charge  density  Ap  so  that 
the  difference  of  these  is  caused  by  the  atomic  core  (nucleus)  or  a  positive  pseudo¬ 
charge  corresponding  to  the  pseudopotential. 

Equation  (4.1)  has  exactly  the  same  form  than  obtained  from  the  pseudopotential 
theory  using  local  pseudopotentials.  The  volume  dependent  term  is  is  the  embedding 
energy  of  an  atom  in  an  homogeneous  electron  gas  without  a  compensating  positive 
background  charge  (subtracting  this  out  gives  the  third  term  in  Eq.  (4.1)).  The  elec¬ 
trostatic  interaction  between  the  atom  and  compensating  positive  charges  enters  now  in 


the  second  sum  as  interaction  between  the  atom  and  bare  pseudopotentials  (which  have 
replaced  the  positive  background).  The  volume  dependent  term  (the  first  sum  in  Eq. 
(4.1))  can  be  cast  to  the  form  familiar  from  the  perturbation  theory  by  noticing  that 
A Ehom  is  the  sum  of  embedding  energies  of  the  bare  pseudopotential  and  free  electrons 
minus  the  binding  energy  of  the  atom  (EATOM).  The  average  electron  density  as 
defined  in  Appendix  C  is  not  exactly  the  overall  average  of  the  valence  electron  den¬ 
sity  which  appears  in  the  formulae  of  the  conventional  pseudopotential  theory.  How¬ 
ever,  if  the  potential  is  weak  the  difference  between  these  two  densities  is  small  and  in 
this  limit  of  weak  pseudopotentials  the  effective  medium  result  agrees  with  the  pseudo¬ 
potential  result. 


Dagens,  Rasolt,  and  Taylor14,15  have  developed  an  method  for  deriving  pair  interac¬ 
tions  in  metals  which  is  related  in  the  effective  medium  theory  in  the  sense  that  the 
results  of  atoms  embedded  in  an  electron  gas  was  used.  In  their  method  the  key  quan¬ 
tity  is  the  induced  electron  density  of  an  atom  embedded  in  a  jellium  vacancy.  This 
has  two  obvious  advantages.  First,  the  electron  density  inside  a  jellium  vacancy 
describes  better  than  the  homogeneous  jellium  the  electron  density  in  the  vacancy  of 
the  actual  metal.  The  second  advantage  is  that  the  density  of  the  electron  gas  is  fixed 
to  the  average  valence  electron  density  which  e.g.  guarantees  that  the  Friedel  oscilla¬ 
tions  have  correct  wave  length.  In  the  theory  of  Dagens  et  al 14,15  the  induced  electron 
density  was  used  to  derive  a  pseudopotential  which  then  was  used  in  normal  way  to 
calculate  the  pair  ion  interaction.  In  the  effective  medium  theory  one  could  also  use 
the  jellium  vacancy  model  as  a  reference  system  instead  of  a  homogeneous  electron 
gas.  One  would  then  first  calculate  the  embedding  energy  of  an  atom  in  a  jellium 
vacancy  and  approximate  the  functional  A £,[«,]  by 


,  ,  ,  ,  Ap(r-R,)Ap(r'-R  ) 

A^KI  =  A  Evactf)  +  - — - J- 

j*i  I  r  ' 


where  the  effective  electron  density  n*'  is  now  defined  requiring 


The  Wigner-Seitz  radius  is  here  determined  from  the  average  density  of  the  metal.  An 
other  possibility  to  satisfy  Eq.  (4.4)  would  be  to  fix  the  density  r%  to  be  the  average 
valence  electn.  n  density  and  vary  i*?^.  In  any  case,  the  use  of  the  systematic  expan¬ 
sion  (2.10)  would  not  any  more  be  straightforward  since  the  evaluation  of  functionals 
A£.[/iy]  would  be  more  difficult.  Also,  the  electron  gas  calculations  would  have  two 
parameters,  the  density  and  the  Wigner-Seitz  radius  which  determines  the  size  of  the 
vacancy.  Approximation  (4.3)  would  nevertheless  lead  to  a  formula  which  would  have 
the  same  form  as  (3.13)  but  where  the  function  F  would  be  related  to  the  embedding 
energy  in  a  jellium  vacancy  (not  necessarily  in  the  same  way  as  in  Fhom  to  A Ehom  in 
Eq.  (3.6)). 


4.2.  Gordon-Kim  model 

Gordon  and  Kim9  have  proposed  a  non  seitcoriMstent  method  for  calculating  binding 
energies  of  molecules.  In  the  simplest  tor:  thi>  theory  makes  an  assumption  that  the 
total  electron  density  of  a  molecule  is  the  si  perposition  of  atomic  electron  densities. 
Density  functional  approximations  are  tl.cn  used  to  estimate  the  energy  change  when 
atoms  are  brought  together  to  form  a  molecule.  The  total  energy  consists  of  kinetic, 
electrostatic,  and  exchange-correlation  pans.  The  electrostatic  contribution  is  exactly 
the  same  as  in  Eq.  (3.13)  when  the  charge  density  Ap  is  taken  to  be  the  charge  density 
of  a  free  atom.  The  simplest  approximation  for  the  kinetic  and  exchange-correlation 
parts  is  a  local  approximation; 

AT  +  A Ex  =  JrfVn(rV(«(r))  -  N\d\na(c)f(na(r)),  (4.5) 


where  n(r)  is  the  total  electron  density  of  the  molecule  of  N  atoms  (supposed  to  be 


similar  to  simplify  notations)  and  /(«)  is  the  energy  per  electron  of  a  homogeneous 
electron  gas  of  density  n.  By  substituting  the  superposition  approximation  (3.1)  to 
n(r),  Eq.  (4.5)  can  be  written  as  follows 


A  T+AEXC  =  XjdVn^r-R,) 


|/(X«a(r-R;))  -f[na( r-R,)) 
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(4.6) 


If  we  define  now  an  effective  density  n?  (different  from  that  used  in  the  normal 


effective  medium  theory,  Eq.(3.9))  as  a  solution  of  the  following  equation 
Jd3ma(r-R,)  /^(r-R.O-^XCr-R,))  -  /(nfl(r-R()-«?*)  |  =  0 

we  can  write  the  total  energy  as 

Etot  =  AJrfV  Ap(r~^)AP(r  ~R/> 

i  ^  i  j*i  'r  ‘  ' 


(4.7) 


(4.8) 


where 


FGK(n)  =  Eatom  +  p3rna(r)^(na(r)+n)-/(na(r))j. 


(4.9) 


Equation  (4.8)  has  now  the  same  form  as  Eq.  (3.13),  but  the  function  F  is  slightly 
different.  Moreover  the  average  density  appearing  in  Eq.  (4.8)  can  not  in  generally  be 
expressed  as  a  simple  superposition  of  (averaged)  atomic  densities  without  doing  addi¬ 
tional  approximations.  The  non-selfconsistent  approach  could  also  be  used  to  estimate 
the  embedding  energy  of  an  atom  in  a  homogeneous  electron  gas.  Since  the  atom 
electron  density  is  kept  frozen  there  will  be  no  Coulomb  contribution  to  the  embedding 
energy  and  it  will  be,  using  the  same  local  approximation  as  in  Eq.  (4.6), 

A Eh™(n)  =  Jci3r[(«a(r)+,-iy(«a(r)+/l)  -  nf{n)  -  ^(r)/(^(r))j.  (4.10) 

The  relation  between  A E^'cn  and  FGK  is  not  exactly  the  same  as  in  Eq.  (3.6)  but  gives 
additional  terms  which  nevertheless  are  smaller  than  the  three  first  term  in  Eq.  (3.6). 
The  use  of  the  local  density  approximation  for  the  exchange-correlation  energy 


functional  in  Eq.  (4.6)  is  commonly  accepted  and  widely  used  in  the  density  functional 
calculations.  For  the  kinetic  energy,  however,  the  local  approximation  (Thomas-Fermi) 
is  not  accurate.  Harris16  has  made  an  extension  to  the  Gordon-Kim  model,  by  calcu¬ 
lating  the  kinetic  energy  change  from  singe-particle  energy  eigenvalues.  This  approach 
clearly  improves  the  local  model,  but  seems  not  to  be  very  useful  here  where  we  want 
to  express  the  energy  with  help  of  the  electron  density  alone.  Recently  Plumer  and 
Stott17  have  made  an  extensive  study  on  approximations  of  the  kinetic  energy  func¬ 
tional.  The  resulting  approximations  for  embedding  energies  of  atoms  in  inhomogene¬ 
ous  electron  gas  could  be  better  used  in  looking  for  approximations  of  the  form  of  Eq. 
(1.2)  for  the  total  energy. 


4.3.  Semiempirical  method  of  Daw  and  Baskes 

Using  the  ideas  of  the  pseudoatom  method  of  Stott  and  Zaremba7  Daw  and  Baskes5 
postulated  that  the  total  energy  is  written  as 

ETot  =  £E(ni)  +  (4.12) 

i  ^  i  j*i 

where  nt  is  the  density  defined  in  Eq.  (3.1)  and  taking  na  to  be  the  free  atom  density. 
The  pair  potential  O  was  assumed  to  be  a  purely  repulsive  short  range  potential 
describing  the  core-core  repulsion.  The  functions  F  and  <t>  where  fitted  to  reproduce 
experimental  lattice  constant,  elastic  constants,  vacancy  formation  and  sublimation 
energies,  and  the  energy  difference  between  FCC  and  BCC  phases.  Equation  (4.12) 
was  found  to  give  good  results  to  surface  geometries  of  solid  Pd  and  Ni5,  and  proper¬ 
ties  of  several  transition  metal  liquids18.  The  comparison  of  Eq.  (4.12)  to  Eq.  (3.13) 
shows  that  even  if  they  have  the  same  form,  the  individual  terms  are  very  different. 
The  pair  potential  in  (3.13)  is  attractive  in  typical  interatomic  distances  in  metals 
whereas  in  Eq.  (4.12)  it  is  weakly  repulsive.  However,  there  is  an  ambiguity  in  deter¬ 
mining  the  functions  F  and  O:  any  part  of  the  function  F  which  is  proportional  to  the 
density  n  can  be  equally  well  included  in  the  pair  potential  sum5,18.  This  means  that 
one  can  always  add  to  the  pair  potential  C>(r)  a  term  y na{r)  where  y  is  any  constant  by 
subtracting  from  the  function  F(n )  the  linear  term  y n.  The  functions  F  and  C>  from  a 
semiempirical  scheme  can  then  not  be  separately  compared  to  those  determined  from 
an  ab  initio  theory. 

In  the  semiempirical  model  one  uses  directly  the  density  without  any  averaging  of 
the  type  of  Eq.  (3.9).  If  one  can  approximate  na~cna,  where  c  is  a  constant,  then  also 
in  Eq.  (3.13)  na  can  be  replaced  by  na  by  redefining  the  function  F.  This  approxima¬ 
tion  is  reasonable  if  only  the  short  range  behavior  of  na  is  important  (meaning  also  that 
one  could  use  the  free  atom  densities  in  Eq.  (3.13)).  The  form  of  the  total  energy 
assumed  by  Daw  and  Baskes  5  can  thus  be  obtained  in  an  approximate  way  also  from 


the  present  theory.  One  should  also  mention  that  in  the  transition  metals  the  semiem- 
pirical  formula  also  includes  the  interactions  between  the  localized  d-electrons,  which 
are  strongly  underestimated  in  the  approximate  formula  (3.13). 

5.  Results 

In  rare  gases  the  function  AEhom(n)  is  proportional  to  the  density  n  and  the  total  energy 
of  the  form  of  Eq.  (3.13)  can  be  expressed  as  sum  of  pair  interactions.  In  Fig.  1  the 
calculated  pair  potentials  for  He,  Ne,  and  Ar  are  compared  to  the  experimental  scatter¬ 
ing  potentials19.  In  the  calculations  the  free  atom  densities  were  used  for  Ap(r).  This 
is  a  good  approximation  in  the  case  of  a  rare  gas  which  is  a  closed  shell  atom  and 
relaxes  only  slightly  when  embedded  in  an  low  density  electron  gas.  The  effective 
medium  theory  can  not  reproduce  the  attractive  part  of  the  potential  which  is  mainly 
due  to  the  van  der  Waals  interaction.  In  the  repulsive  region  the  calculated  potentials 
are  in  fair  agreement  with  the  experimental  results. 

For  A1  metal  we  have  estimated  the  cohesive  energy,  equilibrium  lattice  constant,  bulk 
modulus,  vacancy  formation  energy,  and  surface  energy,  using  for  Ap  both  the  free 
atom  density  and  (in  an  approximate  way)  the  induced  densities  in  homogeneous  elec¬ 
tron  gas.  The  results  are  shown  in  Table  I.  The  free  atom  density  gives  a  too  large  lat¬ 
tice  constant  and  a  too  small  cohesive  energy.  The  reason  is  that  the  function  F(n) 
becomes  too  repulsive  due  to  the  fact  that  the  core  electrons  overemphasized  in  deter¬ 
mining  n.  When  the  atom  is  embedded  in  an  electron  gas  the  screening  makes  the 
atom  more  compact  (only  the  Fricdcl  oscillations  reach  further  out).  This  means  that 
if  Ap  is  the  induced  density  in  an  electron  gas  the  average  density  in  Eq.  (3.13)  will 
be  smaller  and  the  lattice  constant  will  be  reduced  from  the  too  large  value  obtained 
using  the  free  atom  densities.  The  use  of  functions  Ap  calculated  in  a  self-consistent 
way  so  that  it  is  the  induced  density  of  the  atom  when  embedded  in  the  electron  gas  of 


density  «]•  are  very  tedious  to  obtain  due  to  the  long  range  Friedel  oscillations.  As  a 
first  approximation  we  have  completely  neglected  these  oscillations  and  approximate 
that  Ap(r)  is  zero  beyond  a  radius  Rm  which  is  determined  as  the  first  radius  which 
satisfies 

J  J3rAp(r)  =  0.  (5.1) 

r<R„ 

The  radius  Rm  depends  on  the  electron  density  in  which  Ap  is  calculated,  and  for  A1  it 
varies  from  4.5  to  3.9  when  W  goes  from  .001  to  .03.  Since  this  radius  is  smaller  than 
the  interatomic  distance  in  Al,  the  electrostatic  interactions  in  Eq.  (3.13)  reduces  to 
a al(n)n,  where  aat  is  defined  in  Eq.  (B.3)  (the  integrals  go  over  a  sphere  of  radius  Rm ). 
The  total  energy  then  becomes 


Ejot  =  X  2aa/(«i)u,j 


and  depends  in  this  approximation  only  on  the  density.  Since  there  is  no  pair  potential 
part,  the  equation  (5.2)  can  not  reproduce  the  correct  elastic  constants  or  lattice  struc¬ 
ture  (any  structure  would  have  the  same  energy).  However,  it  is  known  that  the  con¬ 
tribution  from  the  pair  potential  sum  to  the  cohesive  energy  for  example  in  Al  is  very 
small3  and  that  the  energy  difference  between  different  lattice  structures,  e.g.  BCC  and 
FCC,  is  small  in  simple  metals.  In  Figure  2  the  total  energy  of  Eq.  (5.1)  is  shown  as  a 
function  of  the  density  n  for  aluminum.  The  cohesive  properties  of  Al  metal  are 
estimated  using  the  expression  (5.2)  and  again  the  free  atom  electron  densities  for 
determining  the  density  n.  The  interpolation  formulae  for  all  the  numerical  data  for  Al 
are  given  in  Appendix  D  and  the  results  are  shown  in  Table  I.  The  agreement  with 
the  experimental  results  is  better  than  in  using  the  free  atom  densities  also  in  as  the 
averaging  function.  This  demonstrates  the  importance  of  the  screening  of  the  atomic 
densities  in  metals.  Still  better  agreement  with  the  experiments  could  perhaps  be 
obtained  by  the  fully  selfconsistcnt  application  of  the  formula  (3.13),  but  the  method 


would  not  any  more  be  practical. 


6.  Discussion  and  conclusions 

Using  the  effective  medium  approach  we  have  shown  that  the  the  total  energy  of  a 
metal  can  be  written  in  the  form  of  Eq.  (1.2),  where  the  function  F  can  be  calculated 
from  the  embedding  energy  of  an  atom  in  a  homogeneous  electron  gas,  and  the  pair 
potential  contribution  becomes  the  electrostatic  interaction  between  screened  atoms. 
The  resulting  equation  (3.13)  reduces  to  the  conventional  pair  potential  result  in  the 
limit  of  a  weak  pseudopotential.  In  the  case  where  the  total  electron  density  is  well 
described  as  a  superposition  of  free  atom  densities,  the  result  of  the  effective  medium 
theory  is  closely  related  to  the  Gordon-Kim  model  of  molecular  binding. 

The  application  of  the  formula  (3.13)  has  some  unwanted  difficulties.  The  charge  den¬ 
sity  Ap  is  the  induced  density  of  an  atom  embedded  in  a  homogeneous  electron  gas  of 
density  n0,  and  has  thus  long  range  Friedel  oscillations.  Since  the  density  7T0  is  gen¬ 
erally  lower  than  the  average  valence  electron  density  of  the  metal,  these  oscillations 
have  a  longer  wave  length  than  the  Friedel  oscillations  in  the  real  metal.  The  existence 
of  these  oscillations  makes  the  summations  in  Eqs.  (3.1)  and  (3.13)  slowly  converging. 
Another  problem  is  the  self-consistency  requirement,  that  n0  should  be  consistent  with 
Ap  through  equations  (3.1),  (3.9)  and  (C.l).  Fortunately,  the  main  part  of  the  electros¬ 
tatic  interaction  comes  from  small  distances  where  the  induced  electron  density  is 
nearly  that  of  a  free  atom  and  fairly  independent  on 

Even  if  the  equation  (3.13)  is  a  result  of  a  systematic  calculation  it  contains  only  the 
first  terms  of  the  exact  expansion  (2.10).  Moreover,  it  has  the  additional  approxima¬ 
tion  (3.1)  of  superposition  of  pseudoatom  densities.  In  the  case  of  closed  shell  atoms, 
rare  gases,  free  atom  densities  can  be  used  to  replace  the  induced  densities  in  the  Eq. 


(3.13).  In  the  case  of  metals,  however,  the  screening  is  more  important  and  can  not  be 
totally  neglected.  Estimations  of  the  properties  of  A1  gives  reasonable  results  for  the 
cohesive  energy  and  equilibrium  lattice  constant  but  underestimates  the  bulk  modulus, 
vacancy  formation  energy,  and  surface  energy. 

In  conclusion,  we  have  derived  an  approximate  formula  for  interatomic  interactions  in 
solids.  In  the  limit  of  weak  pseudopotentials  it  reduces  to  the  conventional  pair- 
potential  formula  for  simple  metals  and  in  the  limit  of  extremely  rigid  atoms  (e.g.  rare 
gas)  it  is  closely  related  to  the  Gordon-Kim  model  for  molecular  binding.  Because  the 
formula  describes  these  two  extreme  limits  correctly,  it  gives  support  for  replacing  pair 
potentials  for  more  general  interatomic  interactions  of  the  form  of  Eq.  (1.2)  in  com¬ 
puter  simulations.  It  seems,  however,  that  ab  initio  calculations  for  the  functions  in 
Eq.  (1.2)  do  not  give  as  accurate  results  as  one  would  want  in  applications  of  the 
method.  A  semiempirical  approach  might  then  be  more  fruitful  in  practical  applica¬ 
tions  to  molecular  dynamics,  and  the  present  theory  would  merely  serve  as  a  guideline 
for  parametrizing  the  functions  needed. 
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Appendix  A 


Assuming  superposition  of  atomic  densities,  Eq.  (3.1),  and  taking  na(Rt-R;)  to  be 
small  if  R,#R;,  we  can  write  Eq.  (3.4)  to  the  second  order  as 

,  ,  0A£(n )  i  d2AE(n.) 

A Ehom(nij)  =  AEhom(rij)  -  na(i-j) — - — —  +  -^(/-y)2  ,  ,  7  (A.l) 


&Ehom(nljk)=  AEhom(nf)  -  [na{i-k)+na(j-k)} 


dd£hom{nk) 

3nk 


1  .  a2 A Ehom(nk) 


+  —  KO'-*)+«a(HO] 


where  a  short  hand  notation  i-y  =  R,-Ry  has  been  used.  Substituting  these  in  Eq.  (2.10) 
keeping  all  term  shown,  we  get 

,  if  dA E(N:) 

Etot=  NEtrost  +  lAE^nj)  -  SMH)  — 

i  2  i  /Vi  ^i 


!  f  )2  a2A£(/If)  1  f  A32A £(«:) 

~  +  -^2  Z'.t;-/)2  — 

6  «  ;Vi  3"«i  12  i  Mi  dnf 


(A.3) 


Since  in  a  solid  each  atom  has  of  the  order  of  10  nearest  neighbors,  it  is  obvious  that 
the  last  term  in  Eq.  (A.3)  is  much  smaller  than  the  other  second  order  term.  If  this  last 
term  is  neglected  we  get  Eq.  (3.5)  by  using  the  definition  of  Eq.  (3.1)  for 

Now  we  demonstrate  that  Eq.  (3.6)  follows  directly  if  we  require  that  the  total  energy 
is  written  in  the  form  of  Eq.  (3.7).  We  assume  a  perfect  lattice  and  derive  the  vacancy 
formation  energy  for  a  unrelaxed  lattice: 

Evac  =  YJLFCriji)  -  F(fii)) 


^  .  dF(ni)  1  ,a2F(/if) 

=^0-0_+_„a0_0_+ 


\w 


wwwm wwwwwwwwwwr^fmwwfrwtwwwwwwwwnwTwrmrrwrTmmm »w%rwu v^- 


-J7  - 
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~n 


dF(ne) 

1  ditj  ' 


(A.4) 


The  second  order  term  can  be  neglected  since  it  is  proportional  to  the  sum  of  squares 
of  individual  atomic  densities  which  is  much  smaller  than  nf.  The  embedding  energy 
is  the  negative  of  the  sum  of  the  vacancy  formation  and  cohesion  energies  and  in  the 
local  approximation  this  means  that 

dF(nt) 


A Ehom(ni)~F(ni)  -  EAT0M  -  nr 


dn t 


(A.5) 


for  which  the  expansion  (3.6)  is  a  solution. 


Appendix  B 

Since  the  electron  density  is  described  as  in  Eq.  (3.1)  the  electrostatic  potential  &f>-<3'(r) 
can  be  written  as 


«OM-r  FnrTr 

j*i  -a.  11  ' 


(B.l) 


We  assume  that  there  is  no  other  positive  charges  inside  a,-  than  the  nucleus  of  the 
atom  in  question.  Then  SV^  is  caused  solely  by  the  homogeneous  background  charge 
of  the  jellium.  In  the  mixed  perturbation  theory  the  effective  density  is  defined  as 


ni  = 


1 


a, 


£  JdVJdV 

at  j*i  a,  cij 


,  Ap(r-R,)Ap(r -Ry) 
ir-r'| 


(B.2) 


where 


f(1  f  ,  Ap(r-R,) 
o„  =  -jj>rf,t>r  |r_r1  ■ 


a  i  a , 


(B.3) 


The  second  term  in  Eq.  (3.14)  can  now  be  written  as 


-tt- 
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_yO 


~aaSlt  =  I  Jd3rJ<fV 


J*i  a,  a, 


,  Apfr-R^Apfr'-R,) 
|r-r't 


(B.4) 


Combining  this  to  the  first  term  of  Eq.  (3.14)  where  the  substitution  (B.l)  has  been 
made  results  Eq.  (3.15). 
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Appendix  C 

In  a  nearly  homogeneous  metal  where  the  pseudopotential  approach  is  good,  the  varia¬ 
tion  in  the  average  densities  «IR,)  at  different  sites  R,  is  small.  We  define 


1  N 

'"o  =  ttITO 
1=1 


and  write  the  first  term  in  Eq.  (3.13)  as 


(C.l) 


XF^RR;))  =  27 


hom 


X«a(R-Ry) 

\j*‘ 


=  s 


Fhom(n0)  +  |X«a(RrRy)  -  "0 


dFhom(n0) 


drio 


(C.2) 


Describing  Fhom  in  terms  of  A Ehom  and  including  only  the  first  two  terms  in  the  expan¬ 
sion  (3.6)  we  get 


£Fhom(F0)  =  I 


_  dAEhom(n0) 


A Ehom(n0)  -  nQ- 


drin 
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1  dAF^ITo) 
~  dnn 


IS 

»  J*‘ 


-ij7^a(RrRr,-) 


,  Ap(r') 


J  lr-r'l 


(C.3) 


where  na  is  described  in  terms  of  na  using  Eq.  (3.9).  Since  a  defined  in  Eq.  (3.10)  can 
be  shown  to  be  exactly7,8 


a  =  + 


9A£^m(n0) 

Sb  ’ 

equation  (C.3)  immediately  gives  Eq.  (4.1). 


(C.4) 


Appendix  D 

The  embedding  energy  A Ehom  was  approximated  by  an  analytic  fit  to  the  numerical 
data.  This  was  chosen  to  be 

A Ehom(n)  =  -1.65  +  1250n  -  490<T120n  (D.l) 

where  n  is  in  atomic  units  and  A E  in  eV.  The  function  aa,(n)  was  approximated  by 

aal(n)  =  -1002  -  68.24/Z2*.  (D.2) 

Both  above  expressions  describe  the  numerical  results  for  aluminum  with  0. 1  eV  accu¬ 
racy  within  the  density  range  «=0.001  •  •  •  0.03.  The  averaging  of  the  atomic  electron 
density  using  the  function  Ap  with  a  cutoff  radius  Rm  defined  in  Section  5  results  to  a 
nearly  constant  enhancement  to  the  local  density.  This  is  approximated  by 

na(r)  =  y in)na{n)  (D.3) 

where  y  depends  on  the  average  density  n  as 

y («)  =  1.369  +  .0478/z1/3.  (D.4) 

The  Gunnarsson-Lundqvist  interpolation  formula-0  for  the  exchange-correlation  poten¬ 
tial  was  used  in  the'calculations  of  A£tom  and  the  electron  densities  of  free  atoms. 


v'VA 


References 


1.  See  e.g.  Interatomic  Potentials  and  Simulation  of  Lattice  Defects  ed.  by  P.C. 
Gehlen,  J.R.  Beeler,  Jr.,  and  R.I.  Jaffee  (Plenum,  New  York  1972). 

2.  W.A.  Harrison,  Pseudopotentials  in  the  Theory  of  Metals  (Benjamin,  New  York 
1966). 

3.  M.  Manninen,  P.  Jena,  R.M.  Nieminen,  and  J.K.  Lee,  Phys.  Rev.  B  18,  7057 
(1981). 

4.  G.  Jacucci,  M.L.  Klein,  and  R.  Taylor,  Phys.  Rev.  B  18,  3782  (1978). 

5.  M.S.  Daw  and  M.I.  Baskes,  Phys.  Rev.  Lett.  50,  1285  (1983)  and  Phys.  Rev.  B 
29,  6443  (1984). 

6.  K.W.  Jacobsen,  J.K.  Nrirskov,  and  M.J.  Puska,  to  be  published. 

7.  M.J.  Stott  and  E.  Zaremba,  Phys.  Rev.  B  22,  1564  (1980). 

8.  J.K.  Ndrskov  and  N.D.  Lang,  Phys.  Rev.  B  21,  2131  (1980). 

9.  R.G.  Gordon  and  Y.S.  Kim,  J.  Chem.  Phys.  56,  3122  (1972). 

10.  H.  Hohenberg  and  W.  Kohn,  Phys.  Rev.  136,  B864  (1964). 

11.  W.  Kohn  and  L.  Sham,  Phys.  Rev.  140,  A 1 133  (1964). 

12.  S.  Lundqvist  and  N.H.  March  (eds.).  Theory  of  Inhomogeneous  Electron  Gas  (Ple¬ 
num,  New  York  1983). 

13.  J.K.  Ndrskcv,  Phys  Rev.  B  26,  2875  (1982). 

14.  M.  Rasolt  and  R.  Taylor,  Phys.  Rev.  B  11,  2717  (1975). 

15.  L.  Dagens,  M.  Rasolt,  and  R.  Taylor,  Phys.  Rev.  B  11,  2726  (1975). 

16.  J.  Harris,  Phys.  Rev.  B  31,  1770  (1985). 

17.  M.L.  Plumcr  and  M.J.  Stott,  J.  Phys.  C  18  4143  (1985). 


18.  S.M.  Foiles,  Phys.  Rev.  B  32,  3409  (1985). 

19.  P.E.  Siska,  J.M.  Parson,  T.P.  Schafer,  and  Y.T.  Lee,  J  Chem.  Phys.  55,  5762 
(1971)  and  J.M.  Parson,  P.E.  Siska,  and  Y.T.  Lee,  J.  Chem.  Phys.  56,  1511 
(1972). 

20.  O.  Gunnarsson  and  B.I.  Lundqvist,  Phys.  Rev.  B  13,  4274  (1976). 


Table  I:  Calculated  properties  of  aluminum  metal.  In  the  non-selfconsistent  model  the 
Eq.  (3.13)  and  in  the  semi-selfconsistent  model  the  Eq.  (5.2)  were  used,  respectively. 
In  both  cases  the  electron  densities  of  free  atoms  were  used.  The  experimental  surface 
energy  is  estimated  from  the  surface  tension  of  the  liquid  phase. 


non-selfconsistent 

semi-selfconsistent 

experimental 

Lattice  constant  (a.u.) 

9.8 

7.3 

7.6 

Cohesion  energy  (eV) 

2.2 

3.3 

3.3 

Bulk  modulus  (Mbar) 

0.1 

0.3 

0.7 

Vacancy  formation  energy  (eV) 

0.05 

0.2 

0.7 

2 

Surface  energy  (111)  (erg/cm  ) 

50 

260 

700 

2 

Surface  energy  (100)  (erg/cm  ) 

30 

170 

700 

Figure  captions 


Figure  1.  Interatomic  potentials  for  He,  Ne,  and  Ar.  The  solid  lines  are  the  results  of 
the  present  calculation  and  the  dashed  lines  the  experimental  results19.  Note  the  change 
of  the  vertical  scale  from  linear  to  logarithmic  at  10  meV. 

Figure  2.  The  variation  of  the  approximate  total  energy  of  Eq.  (5.2)  as  a  function  of 
the  average  electron  density  n  at  a  lattice  site  in  aluminum.  The  quantity  plotted  is 
AF  =  ETOT/N-EATOfl{. 
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